Creation de la base de donnée sur le substance

Dans un premier temps, nous utilisons la fonction « list.files » pour récupérer le nom de tous les fichiers débutant par « bnvd », « bnvd » étant le début commun à toutes les bases de données concernant l’achats de produits phytosanitaires.

Ensuite, nous retirons de cette base de données toutes les lignes où l’identifiant est « 0000 », puis les valeurs manquantes sont transformées en « NA » (ce qui n’était pas le cas avant).

Après, nous enlevons les lignes pour lesquelles la quantité de substance n’est pas renseignée puis, enfin, nous regroupons certains facteurs sous le nom « substance ».

Par la suite, un premier « group_by » est effectué. Pour un même code postal, la même année et la même classification, nous additionnons la quantité de substances achetées.

Le deuxième « group_by » reprend le dataframe obtenu lors du premier « group_by » et permet, cette fois, d’obtenir une seule ligne pour un code postal.

fichiers = list.files(pattern="^bnvd.*\\.csv$")

bddsubstance = do.call(rbind, lapply(fichiers, function(x) read.csv(x, sep=";", dec=".", stringsAsFactors = FALSE,header=TRUE, colClasses=c("code_postal_acheteur"="character"),encoding = "UTF-8" )))

bddsubstance <- bddsubstance[-which(bddsubstance$code_postal_acheteur=="00000"),]

bddsubstance$quantite_substance[bddsubstance$quantite_substance == ""] <- NA

bddsubstance<-bddsubstance[!is.na(bddsubstance$quantite_substance),]

bddsubstance$quantite_substance<-as.numeric(bddsubstance$quantite_substance)
## Warning: NAs introduits lors de la conversion automatique
bddsubstance$classification<-as.factor(bddsubstance$classification)

bddsubstance<-bddsubstance[c(1,2,6,7)]

levels(bddsubstance$classification)[c(4,5)]<-"Toxique"

bddsubstance = bddsubstance %>% group_by(code_postal_acheteur,annee, classification) %>%
  summarise(quantite_totale = sum(quantite_substance), .groups = 'drop')


bddlongsubstance<-bddsubstance %>%
  pivot_wider(
    id_cols = code_postal_acheteur,
    names_from = c(annee, classification),
    values_from = quantite_totale,
    names_glue = "{classification}_quantite-totale-{annee}"
  )

GROUPBY

Voici le résultat du premier « group_by ».

Nous obtenons une quantité totale unique par code postal, par année et par classification.

## # A tibble: 85,752 x 4
##    code_postal_acheteur annee classification quantite_totale
##    <chr>                <int> <fct>                    <dbl>
##  1 01000                 2015 Autre                     44.3
##  2 01000                 2015 N minéral                 30.2
##  3 01000                 2015 N Organique              690. 
##  4 01000                 2015 Toxique                  213. 
##  5 01000                 2016 Autre                    386. 
##  6 01000                 2016 N minéral                479. 
##  7 01000                 2016 N Organique             2671. 
##  8 01000                 2016 Toxique                  274. 
##  9 01000                 2017 Autre                   1460. 
## 10 01000                 2017 N minéral                275. 
## # ... with 85,742 more rows

Le deuxième « group_by » permet de n’avoir qu’une seule ligne par code postal, le transformant en identifiant.

## # A tibble: 6,126 x 17
##    code_postal_acheteur `Autre_quantite-tota~` `N minéral_qua~` `N Organique_q~`
##    <chr>                                 <dbl>            <dbl>            <dbl>
##  1 01000                              44.3                30.2             690. 
##  2 01090                             164.                  7.60           2715. 
##  3 01100                               0.450               2               145. 
##  4 01110                               3.85                4.3             132. 
##  5 01120                             328.                 14.3            8800. 
##  6 01130                               0.00745            NA                24.2
##  7 01140                             392.                 21.4            7377. 
##  8 01150                             314.                 69.4            7123. 
##  9 01160                             115.                 20.1            2391. 
## 10 01170                             396.                  4.76           1119. 
## # ... with 6,116 more rows, and 13 more variables:
## #   `Toxique_quantite-totale-2015` <dbl>, `Autre_quantite-totale-2016` <dbl>,
## #   `N minéral_quantite-totale-2016` <dbl>,
## #   `N Organique_quantite-totale-2016` <dbl>,
## #   `Toxique_quantite-totale-2016` <dbl>, `Autre_quantite-totale-2017` <dbl>,
## #   `N minéral_quantite-totale-2017` <dbl>,
## #   `N Organique_quantite-totale-2017` <dbl>, ...

Création de la base de données sur la quantité de produits phytosanitaires retrouvés dans l’eau

De la même manière que pour la base de données sur les quantités achetées, nous utilisons la fonction « list.files » pour détecter tous les fichiers dont le nom commence par « moy ».

Le dataframe contient les valeurs qui nous intéressent, cependant l’identifiant est « CD_station ». Or, nous souhaitons le code postal.

Pour cela, nous importons « stat », qui nous donne, entre autres, la station (avec son identifiant « CD_station »), son code postal, son code Insee.

Le dataframe « cp », lui, associe à un code Insee le code postal correspondant. Par le biais de deux « merge », nous parvenons à obtenir un dataframe avec les données qui nous intéressent, et le code postal en identifiant.

Enfin, un « group_by » est effectué de manière à avoir le code postal comme identifiant.

fichiers = list.files(pattern="^moy.*\\.csv$")


pe<- do.call(rbind, lapply(fichiers, function(x) read.csv(x, sep=";", dec=",", stringsAsFactors = FALSE,header=TRUE,encoding = "UTF-8" )))
cp<-read.table(file="codepostal.csv",header=TRUE, sep=";", dec=',',colClasses=c("Code_postal"="character"),encoding = "UTF-8")
stat<-read.csv2(file="stations.csv")

stat<-stat[,c(1,2)]
cp<-cp[,c(1,3)]

fus<-merge(pe,stat, by="CD_STATION")

names(cp)<-c("NUM_COM","COD_POS")

bddeau<-merge(fus,cp,by="NUM_COM")
bddeau<-bddeau[,-c(1,2)]




bddlongeau = bddeau %>% group_by(COD_POS,ANNEE) %>%
  summarise(nbpreltot=sum(NBPREL),MOYPTOTAL=mean(MOYPTOT), MAXPTOTAL=max(MAXPTOT), MINMOLRECHTOTAL=min(MINMOLRECH), MINMOLQTOTAL=min(MINMOLQ),MAXMOLQTOTAL=max(MAQMOLQ), .groups = 'drop')

rm(list=c("cp","fus","pe","stat","fichiers"))

Voici les deux dataframes obtenus :

##   ANNEE NBPREL    MOYPTOT MAXPTOT MINMOLRECH MAXMOLRECH MINMOLQ MAQMOLQ COD_POS
## 1  2007      6 0.05333333    0.09        376        377       1       3   01500
## 2  2010      5 0.32900000    0.56        247        409       2       3   01360
## 3  2007      7 0.10028571    0.18        298        377       1       3   01360
## 4  2011      4 0.25000000    0.36        409        409       2       3   01360
## 5  2012      5 0.21720000    0.33        247        409       2       2   01360
## 6  2009      5 0.26260000    0.43        299        383       2       3   01360
## # A tibble: 6 x 8
##   COD_POS ANNEE nbpreltot MOYPTOTAL MAXPTOTAL MINMOLRECHTOTAL MINMOLQTOTAL
##   <chr>   <int>     <int>     <dbl>     <dbl>           <int>        <int>
## 1 01000    2009         8     0          0                  2            0
## 2 01000    2010         4     0.04       0.06             370            1
## 3 01000    2011         8     0.03       0.04             335            0
## 4 01000    2012         8     0.025      0.04             368            1
## 5 01120    2007         2     0          0                  1            0
## 6 01120    2008         4     0          0                382            0
## # ... with 1 more variable: MAXMOLQTOTAL <int>

CARTOGRAPHIE

On transforme le code postal en facteur, puis nous remplaçons « code_postal_acheteur », le nom de la colonne contenant les codes postaux, par « ID ».

Enfin, nous importons le fichier shapefile codes_postaux_région.shp avant d’encoder le fichier en UTF-8.

bddlongsubstance$code_postal_acheteur= as.factor(bddlongsubstance$code_postal_acheteur)
names(bddlongsubstance)[1]= "ID"
codes_postaux= st_read(dsn = "codes_postaux_V5/codes_postaux_region.shp", 
                       layer = "codes_postaux_region",
                       quiet = TRUE) %>%
  select(ID, LIB, DEP)

codes_postaux$LIB= stri_encode(codes_postaux$LIB, from= "ISO-8859-1", to= "utf8")

Nous fusionnons ensuite le fichier shapefile avec la base de données.

fusion= codes_postaux %>% 
  left_join(bddlongsubstance[1:17], by= "ID") %>%
  st_transform(2154)

Enfin, nous créons la carte pour les valeurs « autre » de l’année 2015.

Ce processus est le même pour toutes les autres catégories, il suffit juste de changer le paramètre « col » dans « tm_fill » et « quantité » dans « popup.vars ».

tmap_mode(mode= "view") 
## tmap mode set to interactive viewing
autre_2015= 
  tm_basemap("CartoDB.Voyager")+ tm_shape(shp= fusion)+
  tm_fill(col= "Autre_quantite-totale-2015", palette= "YlOrBr", id= "ID",
          textNA = "Valeur manquante", style = "quantile", n= 6, title= "Quantité totale : <br> Catégorie AUTRE",
          popup.vars = c("Ville"= "LIB", "Quantité"= "Autre_quantite-totale-2015"))+
  tm_borders("black", lwd= 0.3, alpha= 0.6)+
  tm_layout(title = "Quantité de substances phytopharmaceutiques achetées en 2015 (en kilogrammes)", title.position = c("center", "bottom"), legend.bg.color = "white", legend.bg.alpha = 0.4)+
  tm_scale_bar(position = c("left", "bottom"))+
  tm_view(view.legend.position = c("right", "bottom"))

Nous affichons ensuite la carte créée pour la catégorie « autre » en 2015.

TAUX DE CROISSANCE

Pour ce faire, nous calculons, en premier lieu, les coefficients multiplicateurs selon la formule suivante : ((Valeur_n+1 – valeur_n) / valeur_n) + 1.

bdd<-bddlongsubstance

bdd$tauxdecroissancetoxique2016<- ((bdd$`Toxique_quantite-totale-2016`-bdd$`Toxique_quantite-totale-2015`   )/bdd$`Toxique_quantite-totale-2015`)+1
bdd$tauxdecroissancetoxique2017<- ((bdd$`Toxique_quantite-totale-2017`-bdd$`Toxique_quantite-totale-2016`   )/bdd$`Toxique_quantite-totale-2016`)+1
bdd$tauxdecroissancetoxique2018<- ((bdd$`Toxique_quantite-totale-2018`-bdd$`Toxique_quantite-totale-2017`   )/bdd$`Toxique_quantite-totale-2017`)+1


bdd$tauxdecroissanceAutre2016<-   ((bdd$`Autre_quantite-totale-2016`-bdd$`Autre_quantite-totale-2015`)/bdd$`Autre_quantite-totale-2015`)+1
bdd$tauxdecroissanceAutre2017<-   ((bdd$`Autre_quantite-totale-2017`-bdd$`Autre_quantite-totale-2016`)/bdd$`Autre_quantite-totale-2016`)+1
bdd$tauxdecroissanceAutre2018<-   ((bdd$`Autre_quantite-totale-2018`-bdd$`Autre_quantite-totale-2017`)/bdd$`Autre_quantite-totale-2017`)+1

bdd$tauxdecroissanceNorganique2016 <- ((bdd$`N Organique_quantite-totale-2016`-bdd$`N Organique_quantite-totale-2015`)/bdd$`N Organique_quantite-totale-2015`)+1
bdd$tauxdecroissanceNorganique2017 <- ((bdd$`N Organique_quantite-totale-2017`-bdd$`N Organique_quantite-totale-2016`)/bdd$`N Organique_quantite-totale-2016`)+1
bdd$tauxdecroissanceNorganique2018 <- ((bdd$`N Organique_quantite-totale-2018`-bdd$`N Organique_quantite-totale-2017`)/bdd$`N Organique_quantite-totale-2017`)+1

bdd$tauxdecroissanceNmineral2016<-(( bdd$`N minéral_quantite-totale-2016`-bdd$`N minéral_quantite-totale-2015` )/bdd$`N minéral_quantite-totale-2015`)+1
bdd$tauxdecroissanceNmineral2017<-(( bdd$`N minéral_quantite-totale-2017`-bdd$`N minéral_quantite-totale-2016` )/bdd$`N minéral_quantite-totale-2016`)+1
bdd$tauxdecroissanceNmineral2018<-(( bdd$`N minéral_quantite-totale-2018`-bdd$`N minéral_quantite-totale-2017` )/bdd$`N minéral_quantite-totale-2017`)+1

Puis nous calculons le taux toxique moyen annuel pour toutes les catégories, grâce à une moyenne géométrique.

bdd$TAUXtoxiqueMOYENannuel<- ((bdd$tauxdecroissancetoxique2016*bdd$tauxdecroissancetoxique2017*bdd$tauxdecroissancetoxique2018)^(1/rowSums(!is.na(bdd[c(20,21,22)])))-1 )*100   
bdd$TAUXautreMOYENannuel<-((bdd$tauxdecroissanceAutre2016*bdd$tauxdecroissanceAutre2017*bdd$tauxdecroissanceAutre2018)^(1/rowSums(!is.na(bdd[c(23,24,25)])))  -1)*100
bdd$TAUXnorganiqueMOYENannuel<-((bdd$tauxdecroissanceNorganique2016*bdd$tauxdecroissanceNorganique2017*bdd$tauxdecroissanceNorganique2018)^(1/rowSums(!is.na(bdd[c(26,27,28)])))-1)*100
bdd$TAUXNmineralMOYENannuel<-((bdd$tauxdecroissanceNmineral2016*bdd$tauxdecroissanceNmineral2017*bdd$tauxdecroissanceNmineral2016)^(1/rowSums(!is.na(bdd[c(29,30,31)])))-1 )*100   

Cette fonction permet de retirer les valeurs aberrantes de notre dataframe.

bdd <- lapply(bdd[-1], function(x) {
  q1 <- quantile(x, .25,na.rm=TRUE)
  q2 <- quantile(x, .75,na.rm=TRUE)
  IQR <- q2 - q1
  replace(x, x < (q1 - 1.5*IQR) | x > (q2+1.5*IQR), NA)
})

Nous représentons ensuite les taux de croissance moyens entre 2015 et 2018, par catégorie.

bdd<-as.data.frame(bdd)

ggplot(stack(bdd[,seq(29,32)]), aes(x = ind, y = values )) + labs(
  title    = "Taux de croissance annuel moyen de la quantité de produits phytosanitaires achetés ",
  
  x        = "catégorie",
  y        = "taux de croissance annuel moyen entre 2015 et 2018"
) +
  geom_boxplot()
## Warning: Removed 6602 rows containing non-finite values (stat_boxplot).

HYPOTHÈSE FORTE EAU

Dans un premier temps, nous faisons un « group_by » pour n’avoir qu’un seul code postal, puis nous réarrangons l’ordre des colonnes dans le dataframe « longeau ».

longeau<-bddlongeau %>%
  pivot_wider(
    id_cols = COD_POS,
    names_from = c(ANNEE),
    values_from = MOYPTOTAL,
    names_glue = "moyptotal-{ANNEE}"
  )

longeau<-longeau[c(1,6,7,2,3,4,5)]

La même formule que tout à l’heure nous permet de calculer le coefficient multiplicateur.

longeau$TC2008<-((longeau$`moyptotal-2008`-longeau$`moyptotal-2007`)/(longeau$`moyptotal-2007`))+1
longeau$TC2009<-((longeau$`moyptotal-2009`-longeau$`moyptotal-2008`)/(longeau$`moyptotal-2008`))+1
longeau$TC2010<-((longeau$`moyptotal-2010`-longeau$`moyptotal-2009`)/(longeau$`moyptotal-2009`))+1
longeau$TC2011<-((longeau$`moyptotal-2011`-longeau$`moyptotal-2010`)/(longeau$`moyptotal-2010`))+1
longeau$TC2012<-((longeau$`moyptotal-2012`-longeau$`moyptotal-2011`)/(longeau$`moyptotal-2011`))+1

Dans un premier temps, nous supprimons les valeurs non finies et manquantes, puis nous calculons le taux de croissance moyen. Enfin, nous supprimons les valeurs aberrantes.

longeau <- longeau[is.finite(rowSums(longeau[,-1])),]



longeau$TCglobal<-longeau$TC2008*longeau$TC2009*longeau$TC2010*longeau$TC2011*longeau$TC2012


longeau$TCAMmoyenpor<-((longeau$TCglobal^(1/5))   -1)*100



longeauAB <- lapply(longeau[seq(8,14)], function(x) {
  q1 <- quantile(x, .25,na.rm=TRUE)
  q2 <- quantile(x, .75,na.rm=TRUE)
  IQR <- q2 - q1
  replace(x, x < (q1 - 1.5*IQR) | x > (q2+1.5*IQR), NA)
})

Voici le boxplot obtenu pour les taux de croissance de la quantité de produits phytosanitaires retrouvée dans les eaux souterraines :

ggplot(stack(longeauAB[7]), aes(x = ind, y = values )) + labs(
  title    = "Taux de croissance annuel moyen entre 2007 et 2012 en % de la quantité de produits phytosanitaires retrouvée dans l'eau",
  
  x=" taux de croissance moyen en %",
  y        = ""
) +
  geom_boxplot()
## Warning: Removed 37 rows containing non-finite values (stat_boxplot).

LIEN VARIABLE

Nous cherchons à voir s’il y a un lien entre la quantité de produit retrouvée dans les eaux et la quantité de produit acheté dans la première partie nous pour ce faire nous voulons comparer nos variables à un niveau départemental, ainsi, nous ne conservons que les 2 premiers chiffres du code postal, grace à la fonction substr puis nous faisons un group by

bddsub<-bddlongsubstance

eau2012<-longeau[c(1,7)]


bddsub$total2015<-bddsub$`Autre_quantite-totale-2015`+bddsub$`Autre_quantite-totale-2015`+bddsub$`N minéral_quantite-totale-2015`+bdd$N.Organique_quantite.totale.2015
sub2015<-bddsub[c(1,18)]


colnames(sub2015)[1]<-'COD_POS'
cor<-merge(sub2015,eau2012,'COD_POS')


cor$dep<-substr(cor$COD_POS,1,2)

cordep<-cor[c(2,3,4)]
cordep<-cordep[c(3,1,2)]

cordeplong = cordep %>% group_by(dep) %>%
  summarise(Total2015dep=sum(total2015,na.rm = TRUE),moyenne2012dep=mean(`moyptotal-2012`,na.rm = TRUE), .groups = 'drop')


cordeplong<-cordeplong[-3,]
plot(cor$`moyptotal-2012`~cor$total2015,xlab='Total des quantités de produits phytosanitaires achetées en 2015',ylab="quantitée de produit retrouvée dans les eaux en 2012",xlim=c(0,15000),ylim=c(0,0.5)    )

plot(cordeplong$moyenne2012dep~cordeplong$Total2015dep ,xlab='Total des quantités de produits phytosanitaires achetées en 2015 par département',ylab="quantitée de produits retrouvés dans les eaux en 2012",xlim=c(0,300000),ylim=c(0,0.6) )

AGRESTE

Ici nous cherchons à voir s’il y a un lien entre le type de culture et la pollution. Premièrement, nous récupérons la base de donnée de l’agreste, nous enlevons la première colonne qui n’est pas importante puis certaines lignes concernant la corse et les dom tom enfin nous réalisons la matrice de correlation

agreste<-read_excel("agreste.xlsx")


agreste<-agreste[-1]




agreste<-agreste[-c(seq(1,5),100,101),]







cols.num <- c("Céréales","Maïs grain et maïs semence","Oléagineux","Protéagineux","Tournesol","Fourrages et superficies toujours en herbe","Superficie toujours en herbe","Vignes","Jachères")
agreste[cols.num] <- sapply(agreste[cols.num],as.numeric)
## Warning in lapply(X = X, FUN = FUN, ...): NAs introduits lors de la conversion
## automatique

## Warning in lapply(X = X, FUN = FUN, ...): NAs introduits lors de la conversion
## automatique

## Warning in lapply(X = X, FUN = FUN, ...): NAs introduits lors de la conversion
## automatique

## Warning in lapply(X = X, FUN = FUN, ...): NAs introduits lors de la conversion
## automatique

## Warning in lapply(X = X, FUN = FUN, ...): NAs introduits lors de la conversion
## automatique

## Warning in lapply(X = X, FUN = FUN, ...): NAs introduits lors de la conversion
## automatique

## Warning in lapply(X = X, FUN = FUN, ...): NAs introduits lors de la conversion
## automatique

## Warning in lapply(X = X, FUN = FUN, ...): NAs introduits lors de la conversion
## automatique
agreste[is.na(agreste)]<-0


agreste$Département <-substr(agreste$Département,1,2)


names(agreste)[1]<-"dep"

df<-merge(agreste,cordeplong,"dep")



#obligé d'enlever les accents sinon erreur caractère ascii
names(df)<-c("dep","SAU","cereales","ble tendre","Mais","oleagineux","Proteagineux","Tournesol","Colza","Superficie Fourrage","Superficie toujours en herbe","Vignes","Jacheres","Total2015depAchat","Moyenne2012eau")


corPlot(df[-1],upper=FALSE,xlas=3,scale=FALSE,MAR = 1 ,cex=1,cex.axis=1,main="Matrice de corrélation")

Test de moran

au début nous créer les variables contenant la somme des produits phytosanitaires achetés. Au début nous transformons notre fichier en un fichier shapefile puis nous définissons les voisins les plus proches ensuite nous définissons une matrice de poids puis nous pouvons appliquer le test de moran

fusion<-fusion[-c(1,2,3)]

fusion$total2015<-fusion$`Autre_quantite-totale-2015`+fusion$`Toxique_quantite-totale-2015`+fusion$`N minéral_quantite-totale-2015`+fusion$`N Organique_quantite-totale-2015`
fusion$total2016<-fusion$`Autre_quantite-totale-2016`+fusion$`Toxique_quantite-totale-2016`+fusion$`N minéral_quantite-totale-2016`+fusion$`N Organique_quantite-totale-2016`
fusion$total2017<-fusion$`Autre_quantite-totale-2017`+fusion$`Toxique_quantite-totale-2017`+fusion$`N minéral_quantite-totale-2017`+fusion$`N Organique_quantite-totale-2017`
fusion$total2018<-fusion$`Autre_quantite-totale-2018`+fusion$`Toxique_quantite-totale-2018`+fusion$`N minéral_quantite-totale-2018`+fusion$`N Organique_quantite-totale-2018`



fusion_sf<-st_as_sf(fusion)


neighbours_sf <- poly2nb(fusion_sf)

listw <- nb2listw(neighbours_sf,zero.policy = TRUE)


globalMoran <- moran.test(fusion$`Autre_quantite-totale-2015` , listw,na.action=na.omit,zero.policy=TRUE)
globalMoran
## 
##  Moran I test under randomisation
## 
## data:  fusion$`Autre_quantite-totale-2015`  
## weights: listw 
## omitted: 33, 38, 75, 87, 103, 111, 114, 117, 140, 155, 172, 177, 210, 244, 245, 253, 255, 256, 275, 278, 289, 295, 316, 394, 455, 472, 474, 479, 485, 504, 529, 543, 578, 593, 616, 662, 683, 690, 694, 696, 779, 782, 785, 789, 790, 806, 851, 859, 872, 895, 900, 905, 952, 981, 1047, 1056, 1081, 1163, 1165, 1175, 1176, 1341, 1349, 1366, 1368, 1379, 1417, 1430, 1452, 1479, 1489, 1517, 1648, 1666, 1668, 1673, 1679, 1682, 1686, 1692, 1698, 1704, 1705, 1707, 1717, 1718, 1719, 1730, 1731, 1737, 1750, 1751, 1753, 1754, 1755, 1766, 1802, 1803, 1818, 1843, 1846, 1857, 1860, 1862, 1883, 1884, 1892, 1894, 1895, 1897, 1899, 1907, 1929, 1932, 1943, 1952, 1961, 1962, 1965, 2037, 2046, 2061, 2066, 2072, 2078, 2086, 2142, 2156, 2159, 2169, 2171, 2172, 2176, 2177, 2183, 2190, 2192, 2245, 2261, 2273, 2280, 2311, 2314, 2317, 2391, 2447, 2449, 2450, 2457, 2464, 2513, 2521, 2542, 2555, 2564, 2590, 2607, 2611, 2620, 2623, 2625, 2698, 2712, 2728, 2734, 2739, 2811, 2836, 2842, 2854, 2880, 2913, 2921, 2948, 2949, 2950, 2965, 2972, 2983, 2991, 2993, 3003, 3079, 3103, 3144, 3149, 3187, 3197, 3203, 3205, 3207, 3208, 3211, 3237, 3239, 3243, 3244, 3247, 3251, 3254, 3258, 3259, 3271, 3273, 3278, 3284, 3299, 3321, 3336, 3339, 3351, 3357, 3360, 3382, 3389, 3397, 3422, 3424, 3426, 3440, 3441, 3448, 3451, 3461, 3466, 3468, 3469, 3484, 3491, 3524, 3528, 3529, 3532, 3541, 3549, 3551, 3552, 3573, 3578, 3590, 3598, 3606, 3613, 3616, 3622, 3625, 3626, 3636, 3652, 3653, 3695, 3719, 3746, 3747, 3796, 3875, 3885, 3900, 3919, 3932, 3978, 4004, 4016, 4038, 4067, 4068, 4112, 4113, 4121, 4135, 4139, 4174, 4176, 4194, 4199, 4202, 4207, 4208, 4219, 4220, 4221, 4223, 4226, 4232, 4243, 4244, 4263, 4264, 4280, 4327, 4341, 4344, 4348, 4353, 4354, 4355, 4356, 4358, 4361, 4362, 4366, 4369, 4373, 4378, 4379, 4380, 4383, 4384, 4385, 4387, 4393, 4395, 4396, 4397, 4400, 4402, 4404, 4405, 4413, 4416, 4417, 4419, 4428, 4441, 4463, 4466, 4468, 4470, 4483, 4528, 4532, 4537, 4540, 4556, 4557, 4586, 4615, 4643, 4678, 4679, 4707, 4715, 4716, 4720, 4724, 4725, 4726, 4740, 4742, 4745, 4749, 4750, 4751, 4752, 4755, 4756, 4757, 4759, 4764, 4774, 4793, 4837, 4849, 4854, 4860, 4867, 4868, 4869, 4870, 4876, 4921, 4953, 4954, 4961, 4980, 5034, 5036, 5088, 5101, 5118, 5159, 5162, 5164, 5166, 5170, 5200, 5212, 5223, 5225, 5234, 5245, 5272, 5291, 5302, 5311, 5319, 5325, 5327, 5333, 5357, 5380, 5387, 5399, 5406, 5412, 5420, 5436, 5463, 5469, 5474, 5475, 5482, 5497, 5515, 5518, 5537, 5596, 5639, 5642, 5643, 5644, 5646, 5647, 5648, 5650, 5652, 5653, 5654, 5656, 5657, 5658, 5659, 5662, 5664, 5665, 5666, 5669, 5672, 5687, 5689, 5690, 5692, 5696, 5698, 5700, 5703, 5707, 5710, 5714, 5715, 5717, 5718, 5743, 5766, 5774, 5789, 5853, 5862, 5868, 5933, 5957, 5976, 5977, 5980, 6014, 6042, 6043, 6044 n reduced by no-neighbour observations
##   
## 
## Moran I statistic standard deviate = 55.028, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      4.564772e-01     -1.794688e-04      6.886666e-05
local <- localmoran(x = fusion$total2017, listw = nb2listw(neighbours_sf, style = "W",zero.policy = TRUE),na.action=na.exclude,zero.policy = TRUE)

moran.map <- cbind(fusion, local)

tm_shape(moran.map) +
  tm_fill(col = "Ii",
          style = "quantile",
          title = "local moran statistic") 
## Variable(s) "Ii" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.